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Abstract 

A thin plate or slab, prepared so that opposite faces have different surface 
stresses, will bend as a result of the stress difference. We have developed 
a classical molecular dynamics (MD) formulation where (similar in spirit to 
constant-pressure MD) the curvature of the slab enters as an additional dy- 
namical degree of freedom. The equations of motion of the atoms have been 
modified according to a variable metric, and an additional equation of motion 
for the curvature is introduced. We demonstrate the method to Au surfaces, 
both clean and covered with Pb adsorbates, using many-body glue potentials. 
Applications to stepped surfaces, deconstruction and other surface phenom- 
ena are under study. 
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I. INTRODUCTION 



There has been recently a reevaluation of the role of surface stress, as an important 
microscopic indicator of the state of a crystal surface. Recent experimental and theoretical 
aspects of issues related to surface stress have been reviewed by IbachS. 

A good example of connection between micro- and macroscopic quantities is the de- 
termination of surface stress from the measurement of macroscopic deformations of a thin 
plate-shaped sample. According to elasticity theorycl, if a thin plate presents a stress dif- 
ference between its two surfaces, it will bend in order to minimize its total free energy. 
In the limit of small deformations, a simple formula connects stress difference, thickness, 
curvature and elastic properties of the sample. Several methods have been recently used to 
measure the bending. They are mainly optical or based on scanning tunneling microscopy 
(STM)0. Flinn, Gardner and Nixi discussed a laser scanning technique to measure the 
stress-induced curvature, and presented experimental results for stress in Al-Si films as a 
function of temperature. Martinez, Augustyniak and Golovchenkoi measured surface-stress 
changes resulting from monolayer and sub-monolayer coverages of gallium on Si(lll). They 
used LDA calculations of stress in the Ga-covered Si(lll) for the first determination of 
stress in the Si(lll) 7x7 and Si(Ga) superlattice surfaces. Stress variations associated with 
deconstruction of Au(lll) and Au(100) were obtained by STmUB. Moreover, the influence 
of adsorbates on surface stress has been exploited using the sample bending method for 
C/Ni(lll)| S/Ni(lll)i, Co/Ni(100)i, K/Pt(lll)0, Co/Pt (111)0, Ag/Pt(lll)i. 

An atomistic simulation method which controls the curvature of a sample as an extra 
degree of freedom has not been available so far, although this would seem a potentially useful 
tool. The situations in which surface stress varies when physical phenomena occur at inter- 
faces are numerous: adsorbates, reconstructions, steps, phase transitions can change surface 
stress dramatically. All of these phenomena are, under some conditions, well described by 
suitable MD simulations. A simulation method where the slab curvature can be measured as 
a function of adsorbate coverage should represent a powerful connection between important 
microscopic and macroscopic quantities. 

In this work we describe and demonstrate a scheme which is meant to fill this gap. It is 
a classical MD method, but it could in principle be extended to ab initio MD calculations. 
It is based on extending the usual concept of variable cell MD, pionereed by Andersenll! 
and by Parrinello and Rahman0, to a slab with variable curvature. In our formulation the 
curvature is a single, global Lagrangian degree of freedom. Surface stress can be extracted 
as a direct result of the calculation through elasticity equationsEl, which will be described 
below. 

The outline of the work is as follows: section [TJ will present the theory and geometric 
considerations underlying our simulation. In |nj the phenomenology of a bent plate is re- 
viewed in order to extract the pertinent equations. In section [TV] we present some initial 
applications of the method. Finally, section [V| is devoted to discussion and conclusions. 



II. THEORY 

In ordinary variable-cell MD, the coordinates of an atom can be written as 
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where is the position vector of the i-th particle, Sj are scaled coordinates (sf, j = 1 . . .3 
are held between —0.5 and 0.5) and H is a matrix describing the space metric inside the 
cell. Supposing, as was done by Parrinello and RahmanB, that the cell can vary in volume 
and shape, these variations can be accounted for by H, whose elements can be treated as 
an extra set of dynamical degrees of freedom. Our aim is to extend this kind of approach 
to the case of a bending plate, through a different choice of the metric H. 

Let us start with a slab-shaped system, with two surfaces and lateral (x, y) periodic 
boundary conditions (PBC). Suppose we bend the slab cylindrically through a radius R (see 
Figure 1). A convenient choice of the transformation matrix, convenient to our problem, 
could be the following: 

r = H(s,#)-s, (2) 

where H, depending parametrically on the radius of curvature R, is no longer uniform (as 
in the Parrinello- Rahman scheme), but rather depends on the point s. This dependence is 
what originates the curvature. In Figure 1, for example, Si goes from —0.5 to 0.5 along x, 
which we choose to be unaffected by curvature, s 2 goes from —0.5 to 0.5 along direction y 
which follows the curvature, and S3 spans the sample from its inner to its outer surface. 

From now on, we will use as reference surface an ideal "neutral cylinder" inside the slab, 
defined by S3 = 0. In particular, we define the curvature k as the inverse radius of the 
neutral cylinder: 
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(3) 

«3=0 



In this way, if L y is the linear dimension (arc) along the bending direction at k — 0, the 
total bending angle is given by 9m = L y k. Whereas in the Parrinello- Rahman approach 
the whole matrix H is treated as a set of additional degrees of freedom, in our case H is 
completely determined by the single parameter k which we will include in the Lagrangian 
as an extra degree of freedom. 



A. Lagrangian 

The next step is the construction of a Lagrangian. The general, classical Lagrangian for 
a A-particle interacting system is: 

L = T-V = WmAM 2 ~V ( ri ... rjv ) (4) 

1 i 

In Cartesian coordinates ii= H • Sj and ||r|| 2 = s^Gsj, where G = H H is the metric 
tensorS We now define an extended Lagrangian L by including an artificial "curvature 
kinetic energy" : 



L = L+^Wk 2 , 



(5) 



where W acts effectively as a "curvature inertial mass". The velocity of a particle can be 
written as: 

d 



dt v / 



or explicitly 
dt 



(H apS P)=H a p (s,k) sP + Hnpat- 



ds~< ' 



a+ -aT k 



^ + H aP s 



i/3 



(6) 



(7) 



The third term is the usual time derivative for the constant-cell motion, whereas the second 
term contains a time derivative of the extra degree of freedom k. As in the Parrinello- 
Rahman formulationH, we will omit this term from the definition of the particle velocities. 
The kinetic energy of the system will therefore depend on k only through the curvature 
kinetic energy, 



f=-Wk 2 , 
2 



(8) 



whose interpretation is discussed in Appendix A. Parrinello and RahmanEl and AndersenEl, 
using a similar approach, showed that the equations of motion derived from their Lagrangian 
satisfy important requisites, leading in particular to correct thermodynamic averages (not 
dependent on the choice of W), and to a correct balance between external and internal 
stress. 

One can then cast (|7f) as 



OH, 



a P i7 c j8 



OH, 



eta a 



s a + H, 



a/3 



The kinetic energy of a particle then becomes 



-m r 



s a N nf) 8 P = —m s T N s 



-ms iV a/3 



(9) 



(10) 



where N = M -M. 

We need now to express the potential energy in terms of the scaled coordinates and of 
the metric. In the following, we will use mostly potentials (such as pairwise potentials or 
many-body potentials of the 'glue' form) which only depend on pair distances: 

V = V {{Wn - rj 2 / / ./:/../ 1....V}). (11) 

It is therefore necessary to express ||r$ — ij|| 2 in the new coordinates^!. Denoting the metric 
tensor as G (k, Sj) = H T (k, Sj) ■ H (k, Sj) we have 



= r 



'.i 



sjG (k, Si) Si + sjG (k, sj) Sj - 2 sj H T (k, Si ) • H (k, Sj) s r (12) 



We now have all the ingredients needed to write the 3iV + 1 Lagrange equations 

2 m i ~Q~hSi = 



(13) 



k mi N (k, s,) s ; + || - \m x sf^ Sl + mi N (k, s,) s t = 0, I — 1 ... N. 
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B. Construction of the H matrix 



To obtain an explicit form for (^) and ( |T3| ) we write the expression for the Cartesian 
coordinates in a box with linear dimensions (at k = 0) L x , L y and L z : 



x 



s\L x 

V = (1 + s 3 L z ) sin6 

z = ssL z cos 9 + | (cos 9—1) 



(14) 



where 9 is the angle running along the y direction, —0^/2 < 9 < 9m /2 (see Figure 1). With 
the help of simple trigonometric formulas, we can express 9 in terms of s 2 and k: 



■ o ■ kL y 
sin - = 2 s 2 sin — — - 



(15) 



whence equations (|I4]) become: 
- S\L X 



4 s 2 (i + s 3 ^) sin (A; L„/4) ^1 — 4 (s 2 ) 2 sin 2 (k L y / 4) 
s 3 L z (l - 8 (s 2 ) 2 sin 2 (k L y /A)) - { (8 (s 2 ) 2 sin 2 (A; L„/4) 



(16) 



By inspection, we see that a possible choice for H is: 





H(s 2 ,fc) 





V 



o 



I sin(fc Ly/A) 



A / \2 . 2 kL v 

4(s 2 ) sm -J* 



4s2 L,sin^V 1 - 4 ( s 2)^in^ 
L z (l-8( S2 ) 2 sin 2 (^)) ; 

(17) 



ain 2 (fcL y /4) 
" S 2 fc 



2 k Ly 



It can be easily shown that eq. (0) is verified and that, in the limit of zero curvature, 
linifc^o H ( S2, k) = diag [L x , L y , L z ], that is, the Cartesian H matrix is recovered. Moreover, 
the same equations are exact in the limit of finite, fixed curvature. 

The explicit form of N can be obtained by using equations ([H]) and ( |TTD : 



N, 







\ 

q 16(fcL z s 3 +l)^sin 2 (fcn) q 
(l-4(s 2 ) 2 sin 2 (fcti))fe 2 

T2 



(18) 



V Li) 

with u = Ly/A. N is diagonal as a consequence of the chosen cylindrical geometry. 



C. Constrained molecular dynamics 

The Lagrange equations are integrated using numerical methods. In particular, we will 
use a second order velocity Verlet algorithm, which is at the same time well-tested and 
simple. We found some problems due to the instability of the center of mass of the sample 
during the bending, which were overcome and will be commented just below. 
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Our starting geometry is an n— layer crystalline slab with PBC along the x and y di- 
rections. The periodicity is accounted for by the scaled coordinates s% and s 2 . Along s\, 
usual Cartesian PBC are adopted. The boundary conditions for s 2 are also periodic, but 
the images of a point of coordinates (si, s 2 , s 3 ) have coordinates (s%, s 2 ± 1, s 3 ), and lie on 
the cylindrical surface S3 = S3. Every point of the sample is equivalent along s 2 and so a 
"translational" invariance is still present, albeit in a curved geometry. 

For k 7^ the dynamics generated by ([HJ) does not conserve the position of the center of 
mass of the slab along the radial direction S3. To keep the center of mass fixed, we resort to 
a constrained MD scheme. Constrained MD is a well known technique, and we adopted the 
method used by Ryckaert, Ciccotti and BerendserJUl and Andersen@ for simulations of rigid 
molecules. The main point to be kept in mind is that constraints must be verified exactly 
at every time step, otherwise instabilities will occur. Technical details on the practical 
implementation of this constraint are given in Appendix B. 

The general MD protocol we apply consists of either constant energy runs or constant 
temperature runs (up to 300 K) obtained by velocity rescaling, followed by quenching to 
T = 0. Occasionally, we also introduced a damping term in the curvature dynamics. We 
chose a typical time step of the order of 10~ 14 s. 



III. PHENOMENOLOGY OF A BENT PLATE 

Before moving on to the actual implementation, it is useful to recall some known results 
concerning the phenomenology, largely contained in the review by IbachS. 

The definition of surface stress, according to GibbsS, as the "reversible work per unit area 
required to stretch a surface elastically" , points out the difference of this physical quantity 
with respect to the surface free energy, defined as the "reversible work per unit area to 
create a surface". Stretching a solid surface implies modifying the substrate, whereas a 
simple increase of area does not: whence the difference. 

Let us start with the bulk stress tensor, whose element Ty is defined as the i-th component 
of the force per unit area acting on the side (with the normal to the surface parallel to the 
j'-th direction) of a small cube in the sample. The corresponding surface stress tensor can 
be defined as the deviation of this quantity with respect to the bulk value, integrated along 
the surface normal: 



Ta{z)-T$ (19) 



The integral gives non-zero contribution only where the value of the stress deviates from 
the bulk value, i.e. at the surface. 

Let us consider a rectangular shaped sheet of thickness t, delimited by two identical (100) 
surfaces. For our purposes the following simplified construction can be put to use without 
loss of generality!. Schematize the plate at zero curvature as a "substrate" of thickness t and 
length L along X\, plus two films adsorbed on it: film A on the lower part of the plate, and 
film B on the upper part. Let La the unconstrained length (along x±) of film A, and Lb the 
length of film B. Let Lb > L and La < L. In order to match the length L of the substrate, 
the film must be stretched (compressed) by an amount ALj = L — Lj, with I = A,B. 
The deformed film attached to an undeformed substrate does not represent a condition of 
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minimum free energy; the free energy of the entire system can be lowered by deforming the 
substrate slightly so as to reduce the deformation of the film. This deformation has two 
components: an overall compression of the substrate and a bending of the substrate. 

It can be shown! that the neutral plane (i.e. the plane along which the strain is zero) can 
be set approximately in the middle of the sheet, and the deviations from this approximation 
are negligible for most cases. In the following, we will extract the conditions for the minimum 
free energy of the sheet, allowed to bend along a single direction. For simplicity, we will 
choose this to be the [100] direction (of course similar formulae can be derived for any other 
choice) . If we denote the coordinates along the bending direction, the non-bending direction 
and the normal to the surface with si, S2 and S3 respectively, the situation is similar to the 
one depicted in Figure 1. Each element of the sample is subject to a strain en (£3). The 
change in free energy per unit area associated to this strain, of an element at height S3 is: 

reii( S3) 

u{s 3 )= / r n (s 3 ) den. (20) 







The total free energy change of the bent sample (per unit area) is: 



U = I ds 3 / r n (s 3 ) de n . (21) 



2 



We can separate this integral into two surface parts and a bulk part: 

U = U {s+) + U (s - ) + U {b) (22a) 
= / +l ds 3 



l/<->= ds 



'+/3 feu («s) 



ruC^-rff de al (22b) 



! 3 
2 







I \ ( 6 ) 
Til ( S3) - 7lY 



den (22c) 



= [** ds 3 f ll(S3) rff ( s 3 ) de n . (22d) 



_> 







The first two terms lead to the definition of surface stress in the Gibbs sense, and can be 
written with good approximation: 



U^ + U^ = (t$(+)-tI?(-)) ^ , (23) 



where ±kt/2 are the strains at the upper and lower surfaces. 

The bulk term can be expressed using the elastic constants of the crystal. Hooke's law 
relates strain and bulk stress tensor components in the following way: 

£ll = (Til Til + 012 T 22 ( 24) 

£22 = 012 r ll + 011 r 22 • 

In our case, e22 = 0, therefore, after introducing Young's modulus Y and Poisson's ratio v: 

Y=±- (25) 
011 



v - 



012 
011 
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one obtainsS: 



£11 = y ^"ii- (26) 



If we substitute fl26|) into ( |22d|) we obtain, after integration, 



um = iiThs- k2t3 ' (27) 



the total free energy of the bent plate is then: 

M^M-^H^yJ+i^^ 3 . (28) 

Minimizing F with respect to the curvature, one obtains the equation relating the equilib- 
rium curvature hu to the difference in surface stress between upper and lower surface: 

^""SRib)*"* 1 (29) 

where Y is the Young's modulus, v the Poisson's number, kyi the curvature and t the 
thickness of the sample. This relation is known as Stoney's equation, and was first derived 
in 1909 (except for the biaxial nature of the stress@'0). We stress that this formula is valid 
for a uniaxial bending, whereas for an isotropic bending in two directions the factor (1 — u 2 ) 
should be replaced by (1 — v). Another limit of validity is that the plate must be infinitely 
large, relative to its thickness. In other words, there will be finite thickness corrections 
in (t/L) where L is the lateral (x,y) size. We did not make any attempt to find these 
corrections so far. 

In the remainder of this section we present the simple derivation of the elastic constants 
two high symmetry surfaces: (100) bent along [Oil], and (111). 

We start from the elastic stiffnesses CV,- which are related to the compliances by the 
equations^: 

a ^ Cii + C12 

(Cu — C12) (Cu + 2 C12) 



C 



C12 
044 



12 



(C11 — C12) (C11 + 2 Cu) 

1 

C44 



If the uniaxial bending is not along a [100] direction and for surface orientations other than 
(100), the compliances an and a 12 which appear in ( p5|) have to be replaced by effective 
elastic constants, which in the (100) case, with bending along a [011] direction, have the 
form: 



1 / 1 



a 
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0"l2 = °"12 + £ \ °"H + °"12 + 7^°"44 



o"ii - - ( °"n + °"i2 + 2 a44 ) ( 31 ) 

1 / 1 



s 



whereas in the (111) case have the form: 

1 



a 
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2 V 2 



o"n - - ctu + (Ti2 + -a 44 (32) 



1 / 1 

^ 12 = Cl2 + - I (Tn + Cl2 + - a 44 



IV. APPLICATION TO METAL SURFACES 
A. Implementation 

The first goal of the variable-curvature MD simulation is to obtain an equilibrium value 
for the curvature in order to extract the surface stress difference, according to equation (p9|) . 
We can then compare the outcome with the surface stress difference calculated independently 
using Kirwood's formula!!, and assess the success of the method. We will also verify the 
validity of the approach by analyzing the behavior of the curvature as a function of the slab 
thickness. 

We will present two different exemplifications. Far from being exhaustive, these results 
are only meant to show the feasibility of the method in view of later applications to surface 
physics. 

First we simulated several Au slabs with (100) orientation and different types of recon- 
struction on the two sides. This test can show the sensitivity of the method: differences 
of the order of 5 me v/A 2 can be appreciated. Next we examined the effect of isotropic or 
anisotropic adsorbates on the surface stress, as exemplified by Pb/Au(100). 

All simulations were performed using many-body potentials of the glue typeil. They 
have the general form: 

V = 5E*(r<,)+£l7(rO (33) 

i, j i 

j 

where rii is a generalized atomic coordination. $(r), U (n) and p (r) are empirically 
constructed, by fitting several properties of the system. Well tested glue potentials are 
available!! for Aui and Pbi. This scheme can be extended to binary systems by intro- 
ducing a mixed two-body potential $ab(?"), and assuming rii to be a linear superposition of 
density contributions pa and ps supplied respectively by A- and I?-type atoms to site 



v= l - 

2 



E E *AA (%) + E *AB (ry) + E E ®AB (r y ) + E ®BB 
ieA \jeA jeB I ieB \jeA jeB 



+ Y,Ua {ni) + Y.U B (ni), (34) 

i£A ieB 



with 



Hi = E W aPa (r y ) + E W bPb (r y ) . (35) 
j£A jeB 

where Wa and Wb are suitable weights. In our case W\ u = 1.045 and Wpt, = 0.957. 
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B. A clean surface: change of surface stress with reconstruction 

We prepared a family of Au(100) slabs composed by about 600 atoms per layer (L x = 
28.78 A, L y = 172.69 A), and different thicknesses of 8, 12, 16 and 20 layers, with [Oil] as 
the bending direction. 

In order to extract numerical values for the surface stress from equation fl29"|), we need 
the correct elastic constants for this case. For the case of Au (as described by the glue 
model), the elastic stiffnesses have been calculated! 2 ^, giving 

C u = 2.203- 10 7 N/cm 2 (36) 
C12 = 1.603- 10 7 N/cm 2 
C 44 = 0.600- 10 7 N/cm 2 . 

With these values, using (p5|) and (0), the Young's modulus and the Poisson's ratio for 
Au(100) , with bending along [011], are F (lo0 ) = 1.322 ■ 10 7 N/cm 2 and z/ (100) = 0.1022. 

Both (100) surfaces were prepared in the reconstructed state, characterized by a close- 
packed triangular overlayer on a square substrate. We chose a fixed reconstructed (1x5) 
structure (six [011] rows on top of five) for the reference (lower) surface of the plate, and 
we changed the nature of the upper surface, ranging from (1x3) (i.e. four surface rows 
over three substrate rows, a "compressive" situation) to a (1 x 20) structure, (a "tensile" 
situation), touching the (1 x 4), (1 x 5), (1 x 6), (1 x 10), (1 x 12), and (1 x 15) structures. 
In the (1x5) case no bending was expected or observed, since the surfaces on the two sides 
were identical. We also considered (11 x 3), (11 x 4), (11 x 5),... obtained by adding [011] 
rows. The lowest energy surface in the glue model is (34 x 5)e3, in the experiment is close 
to (28 x 5), but all (M x 5) surfaces differ very slightly in free energyn3lH3. The essential 
point here will be that surfaces which differ only very little in energy, for example (1 x 5) 
(102.3 meV/A ) and (1 x 12) (103.2 meV/A ) may differ enormously in surface stress. 

A typical snapshot of the simulation is shown in Figure 2. 

Figure 3 shows the time evolution of the curvature k. Starting from zero, the slab 
curvature reaches the equilibrium value very rapidly. A damping term, although not really 
necessary, has been added in order to speed up the convergence. 

We verified that, as it should be and as can be seen from the figure, the value of the 
mass W is irrelevant in determining the equilibrium value of k. In the (1 x 3) and (1 x 
4) cases, the equilibrium k was positive, whereas in the other cases it was negative. The 
simulation was performed at low temperature (from a few degrees to 300 K). The difference 
in surface stress between upper and lower surface has been extracted using (^9[). The surface 
stress of the two isolated surfaces has been independently obtained from MD forces using 
the standard Kirkwood-Buff formulall. The comparison of the stress differences is excellent, 
as can be seen from Figure 4, showing that the variable curvature method works. The only 
deviation is for (1x3), and clearly attributable to excessive curvature. 

The calculation is repeated with a variety of different unit cells in the top surface; by 
adding a [011] row every ten, we obtain (11x3), (11x4), (11x5),... unit cells instead 
of (1 x 3), (1 x 4), (1 x 5), ... The lower (1 x 5) surface of the slab was left unchanged. 
The stress variations are cleary very large, both relative to the basic stress of (1 x 5) 
(203.12 meV/A 2 ) and to the surface energy (~ 100meV/A 2 ). We judge that the smallest 
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surface stress difference detectable with our kind of simulation is of the order of 10 me v/A 2 . 
This sensitivity should be very important in the study of phase transitions. 

We have also verified the linear dependence of the curvature from the inverse square of 
the thickness predicted by equation (^). Figure 5 shows the fit with 8, 12, 16 and 20 layers, 
and the agreement is rather good. 

C. Change of surface stress with adsorption 

The surface stress is obviously a strong function of adsorption. As an example, we 
will consider here Pb on Au. Underpotential Deposition (UPD)i has been recently used 
together with STM techniques to exploit the formation of 2D phases of Pb on a substrate of 
Ag(100) and Au(100). Through different electrochemical potentials, one can obtain different 
phases, such as Au(100)-c(2 x 2) Pb, Au(100) c (3^ x y^) R45° Pb, Au(100) c(6 x 2) Pb. 
In particular, a close packed layer of lead shows a slight contraction in both [Oil] directions 
of the quadratic substrate lattice; the final structure is pinned to the substrate lattice at 
equivalent adsorption sites, leading to an Au(100)-c(6 x 2) moire superstructure. 

Using the potential in eq. (|34|), we found precisely this c(6 x 2) superstructure as a 
local energy minimum, through quenching of a sample prepared with a perfect substrate of 
Au(100), and a triangular overlayer of Pb with the same density of a bulk (111) layer. In 
the optimized structure, the (6 x 2) pattern is evident (Figure 6(a)). Domain walls appear 
as the result of overlayer contraction. We were able to obtain a perfect (6 x 2) (Figure 6(b)) 
by starting from an initial configuration with a denser overlayer. Deeper energy minima are 
obtained by the simple commensurate square overlayers with vacancies as in the case shown 
in Figure 6(c). 

We performed variable curvature simulations on both the imperfect and the perfect 
c(6 x 2), and on the square overlayer with vacancies, with Au(100) (1 x 4) as the reference 
(t 22 = 147.52 meV/A 2 ). For c(6 x 2), which is anisotropic, two different simulations were 
carried out in order to access separately r n and r 22 (the adlayer was rotated of 90 degrees, 
so as to keep r 2 2 as reference). 

For the (6 x 2) structure with domain walls we found r 22 = 88.10 me v/A 2 {t u was 
not measured in this case), and for the perfect (6 x 2) structure 7~n = 72.41 meV/A 2 and 
r 22 = 38.52 me V/A 2 . For the square structure with vacancies, we found r 22 = 40.34 meV/A . 
The stress anisotropy is therefore large, and the effect of the domain walls is also remarkable. 

V. CONCLUSIONS 

In this work we have introduced a variable curvature slab simulation method which 
represents a promising technique to study properties of crystal surfaces, in particular surface 
stress differences. 

The first tests of the method, based on the comparison with known results, reveal a good 
accuracy, in the range of 10 meV/A , and easy applicability. We are now considering studies 
of a variety of phenomena, including surface phase transitions, using this approach. 
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APPENDIX A. INTERPRETATION OF THE KINETIC ENERGY OF THE 

CURVATURE. 

If the arc S3 = is forced to have length L y during the whole simulation, the volume Q 
of the sample is guaranteed to remain constant and equal to L x L y L z . This is easily seen by 
performing the integral in polar coordinates. 

The difference in surface area between upper and lower layer is of course: 



AA 



MAX 



R + 







MAX 



R 



2 



The relative time variation of this area is 

d /AA\ 



V d 



dt \ A J L x L y dt 



Therefore, the fictitious kinetic energy term is: 



x^z u MAX 



k — L ? k 



L x L y L z k 



fife. (37) 



(38) 



T 



-Wk 2 





(AA) 



(39) 



This result is simple but important because allows a connection with Andersen molecular 
dynamicst3: whereas in Andersen's formulation the extra degree of freedom is the volume 
Q and the term ~W is the kinetic energy tied to volume variations, in this case there is a 
variation in the area difference between the two faces of the sample. The analogy could be 
brought further, adding to our Lagrangian a term —Pk, which is the reversible work of a 
radial force decaying as the square of the curvature radius, acting identically on all particles, 
whereas the term — P Q in Andersen's formulation is the reversible work of the hydrostatic 
pressure during the variation of the cell volume. Such an external force, acting as a torque, 
could be used to force the equilibrium curvature to a different value from the one suggested 
by stress imbalance between the two surfaces. Applications of this point are presently under 
study. 



APPENDIX B. CONSTRAINED DYNAMICS FOR FIXING THE CENTER OF 

MASS 

A. The velocity Verlet as a predictor-corrector. 

Following^, it is possible to cast a particular velocity Verlet algorithm as a predictor- 
corrector algorithm. 
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One can start from 



r (t + St) = r (t) + rj (t) + r 2 (t) (40) 
Tl (t + St) = n (t) + r 2 (t) + r 2 (t + St) 

where, as usual, 

1 H n r 

r -W = SfW"dF- (41) 

This can be written as a two stages predictor-corrector, with a force evaluation in between: 

r (t + St/2) = r (t) + n (t) + r 2 (t) (42) 
r 1 (t + St/2) = r 1 (t)+r 2 (t) 

r (t + St) = r (t + St/2) 

n(t + <ft) = n(t + 5t/2) + r 2 (t + 5t) 

where r 2 (t + 5t) is obtained by the equations of motion. In this case the predicted values for 
the positions coincide with the corrected ones, whereas the velocities are corrected (second- 
order, or three values, predictor-corrector). 

When constraints are present, the forces acting on the i-th particle can be written as: 

Fi = 5, + g,«t + g! J) . (43) 

The constraints forces g a at each stage are approximated with gjf\ where (x) = r,v, so as 
to ensure that the positions and the velocities, respectively, satisfy the constraints, and the 
quantities at the various stages have to be corrected to keep into account the constraints 
forces. 



B. The bent plate case. 

In the case of varying curvature molecular dynamics, the bending leads to a drift towards 
the center of curvature of the center of mass. In order to fix the S3 component of the center 
of mass in the origin, the only constraint to be satisfied is (assuming all the particles to have 
the same mass for sake of simplicity) 

N 

a = 5> 3 , = (44) 

i=l 

where are the scaled coordinates in our curvilinear system (in this appendix, s, s and s 
are normalized according to ([IT])). The time derivative of the constraint equation gives a 
constraint on velocities: 

1 N 

*=7lE** = 0. (45) 

N N 

The Lagrangian must then be completed with a term A# J2 s 3i + V E &3i- 

i=l i=l 
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Two different approximations to the constraint forces have to be chosen so that positions 
and velocities satisfy the constraints exactly. In this particular case the forces turn out to 
be simple constants: 



$ = -A« (46) 

(v) 

The first equation, involving positions, is: 



n (v) - A 
<?3i — ~ A V 



S 3i (t + St) = S 3i (t) + S 3l (t) + S 3 i (t) - ^~^R (47) 

Xr has to be chosen in order to satisfy the position constraint at time t + St.: 



££i (s 3i (t) + s 3i (t) + s 3i (t) - Xr) = 

x _ r> Efal(*3i(*)+i3i(*)+g3i(*)) 

Xr - Zrm . 



Therefore the former equation becomes: 



(48) 



s 3 i{t + St) = s 3i (t) + s 3i (t) + s 3i (t) J - — . (49) 

The second equation, involving velocities, is: 

s 3i (t + St) = s 3 i (t) + hi (t) - l^-^R + hi (t + St) - V (50) 



(51) 



Xy must satisfy 

Ei=i I S 3l (t) + S 3i (0 n + S 3* (* + St) ~ ^-Ay I = 

Ay - Zrm . 

and the equation becomes 

i (f + XA i (A + S (A ^=l( g 3. (t) + S 3i (t)+S3i (*)) 

s 3i (t + St) = s 3i (t) + s 3i (t) — (52) 

^ (.^m J2f =1 (hAt + St)-s 3j (t)) 
+s 3i {t + bt) J - — . 
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FIGURE CAPTIONS 



Figure 1. The geometry of a bent plate with curvilinear coordinates; s x and s 2 lie 
in the neutral plane (shown as lighter in the figure) whereas s 3 is vertical. R is the 
curvature radius measured at s 3 = 0, and 6 M is the bending angle. 

Figure 2. A typical snapshot of the simulation, with a positive value of the curva- 
ture. The top surface is a (1 x 10) and the lower surface a (1x5), both hexagonally 
reconstructed Au(100). The corrugations of both surfaces are evident. This particular 
sample will invert its curvature during the simulation, and the equilibrium value of k 
will become negative. 

Figure 3. Time evolution of the curvature k for two different values of the 'mass' W. 
On the right axis, the bending angle 6m is shown. 

Figure 4. Comparison between variable curvature simulations and Kirkwood-Buff 
formula, with different structures at the upper surface of the sample. The lower surface 
is a (1x5) in every case. Right ordinate axis: equilibrium value for the curvature. 
Left ordinate axis: corresponding difference in surface stress, according to Stoney's 
equation. In order to obtain the absolute value of the surface stress component r 2 2 of 
the upper surface, the reference (1 x 5) result (203.12 meV/A 2 ) has to be added. 

Figure 5. Dependence of the curvature upon thickness in the case of (1 x 10)/ (1 x 5). 
The inverse square dependence predicted by Stoney's equation is recovered, and the 
stress difference is extracted from the fit. 

Figure 6. Different structures for the Pb overlayer (grey and black atoms) over a 
square Au(100) substrate (white atoms), (a): (6x2) structure with domain walls. 
Black atoms are higher, two MD boxes are shown along x; (b): (6 x 2) without domain 
walls; (c) square adlayer with vacancies. 
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